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Coarsening  of  the  nickel  phase  is  known  to  occur  in  solid  oxide  fuel  cell  (SOFC)  anodes  consisting  of  Ni 
and  yttria-stabilized  zirconia  (YSZ).  However,  the  exact  nature  of  the  coarsening  process  is  not  known, 
nor  how  it  affects  three-phase  boundaries  (TPBs)  and  the  resulting  electrochemical  performance.  We 
apply  a  phase-field  approach  to  simulate  the  microstructural  evolution  of  Ni-YSZ  anode  functional 
layers.  An  experimentally  obtained  three-dimensional  reconstruction  of  a  functional  layer  from  an  anode- 
supported  SOFC  is  used  as  the  initial  microstructure.  The  evolution  of  the  microstructure  is  characterized 
quantitatively  by  examining  the  TPB  density,  interfacial  area  per  unit  volume,  and  tortuosity  versus  time. 
The  assumed  TPB  contact  angles  are  found  to  have  a  strong  effect  on  the  microstructural  evolution;  in 
particular,  reducing  the  contact  angle  of  nickel  on  YSZ  yields  less  TPB  reduction. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Nickel  particle  agglomeration  and  coarsening  have  been  con¬ 
sidered  major  factors  responsible  for  degradation  in  contemporary 
SOFCs  with  nickel-yttria-stabilized  zirconia  (Ni-YSZ)  anodes  [1]. 
Despite  its  importance,  coarsening  is  difficult  to  investigate  sys¬ 
tematically  because  it  proceeds  slowly  and  must  be  examined 
over  a  long  period  of  time.  A  few  experiments  of  Ni  coarsening 
can  be  found  in  the  literature  [2-4],  but  they  do  not  examine  the 
microstructures  in  three  dimensions  (3D),  and  thus  it  is  difficult  to 
characterize  the  evolution  of  tortuosities  or  connectivities  of  vari¬ 
ous  phases  [5].  There  have  been  some  modeling  efforts  as  well,  but 
these  models  are  either  statistical  [6]  or  semi-empirical  [3,7,8]. 

Ni  coarsening  in  the  SOFC  anode  is  a  capillarity-driven  phe¬ 
nomenon.  Regions  with  high  curvatures  have  higher  chemical 
potentials  than  those  with  lower  curvatures  in  accordance  with  the 
Gibbs-Thomson  effect,  and  thus  the  material  will  be  transported 
from  these  regions  to  lower-curvature  regions  when  mobility  is 
sufficiently  large  for  the  time  scale  of  interest.  This  leads  to  lower 
free  energy  of  the  system.  The  process  can  be  formulated  as  a  sharp- 
interface  free-boundary  problem,  but  it  is  difficult  to  numerically 
solve  in  3D  since  explicit  tracking  of  the  interface  is  required.  The 
phase-field  model  (PFM)  handles  the  interface  implicitly,  where  the 
interface  information  is  embedded  in  field  variables,  often  referred 
to  as  the  order  parameters  (OPs).  Immiscibility  and  interfacial 
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energy  are  naturally  incorporated  in  the  PFM  through  the  bulk  free 
energy  and  the  gradient  energy  penalty  over  the  diffuse-interfacial 
region.  The  interfacial  energy  ratios  among  different  interfaces  can 
be  specified  in  the  PFM  through  proper  parameterization  of  the  free 
energy  functional. 

In  this  paper,  we  present  two  models  to  quantitatively  sim¬ 
ulate  Ni  coarsening.  Our  models  are  free  from  empirical  fitting 
parameters  when  input  kinetic  and  thermodynamic  coefficients  are 
known.  Furthermore,  to  the  best  of  our  knowledge,  this  work  is  the 
first  to  perform  coarsening  simulations  in  3D  using  an  experimental 
microstructure  of  a  Ni-YSZ  anode. 


2.  Model  formulation 

To  simplify  our  models,  a  few  approximations  have  been  made. 
First,  the  volume  of  each  phase  in  the  Ni-YSZ  anode  is  assumed  to 
be  conserved  (which  is  equivalent  to  mass  conservation  since  the 
material  is  incompressible).  (Note  at  ordinary  operating  tempera¬ 
tures  of  a  SOFC  (500-1000  °C),  the  evaporation  of  Ni  is  negligible 
[9].)  Therefore,  we  apply  the  Cahn-Hilliard  (conserved)  equation 
as  the  governing  equation.  Second,  the  kinetic  and  thermody¬ 
namic  properties  are  estimated  at  1 1 00 0  C  at  which  complementary 
experiments  are  conducted.  In  addition,  Ni  surface  tension  is 
assumed  to  be  independent  of  crystal  orientation,  because  its 
anisotropy  is  insignificant  at  this  temperature  [10].  For  the  surface 
diffusivity  of  Ni,  we  adopt  the  value  measured  from  the  mass  trans¬ 
fer  method,  which  can  provide  the  mean  value  of  multiple  crystal 
facets  and/or  the  value  of  certain  facet  of  a  single  crystal  [11]. 
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We  further  simplify  our  models  by  identifying  the  dominant 
transport  mechanism  that  dictates  Ni  coarsening  in  a  Ni-YSZ  anode. 
Three  mechanisms  play  a  role  in  the  evolution  of  Ni  microstructure: 
evaporation-condensation,  bulk  diffusion,  and  surface  diffusion. 
The  evaporation-condensation  mechanism  is  negligible  here  due  to 
the  very  low  Ni  vapor  pressure  [9].  Assuming  atoms  on  the  interface 
have  similar  values  of  chemical  driving  force  for  both  surface  and 
bulk  diffusion  (a  quasi  surface-volume  equilibrium  assumption), 
the  ratio  between  surface  diffusion  and  bulk  diffusion  contributing 
to  the  transport  can  be  represented  as  o  =  Ds8s/DBw,  where  Ds  and 
Db  are  surface  and  bulk  diffusivities,  respectively.  Surface  diffusion 
is  considered  to  occur  within  the  first  several  lattice  layers  on  the 
surface  within  a  thickness  8S  ~  1  nm.  In  the  anode  microstructural 
data  we  acquired,  the  characteristic  length  w  of  Ni  particles  falls 
between  0.1  and  1  p,m.  According  to  the  literature,  Ni  surface  diffu- 
sivity  is  on  the  order  of  10-10  m2  s-1  [10,1 1  ],  while  bulk  diffusivity 
is  of  10-16  m2  s-1  at  1100  °C  [12].  Therefore,  the  surface  to  bulk  dif¬ 
fusion  ratio  is  cr^  103  in  this  system,  which  suggests  that  surface 
diffusion  dominates. 

It  is  generally  believed  that  the  YSZ  phase  has  very  low  mobil¬ 
ity;  that  is,  YSZ  serves  as  the  supporting  structure  in  which  Ni 
coarsens.  Therefore,  we  may  either  allow  YSZ  to  transport  with  its 
low  mobility  (Model  A),  or  assume  it  to  be  immobile  (Model  B). 
The  framework  of  these  two  PFMs  will  be  discussed,  and  we  will 
demonstrate  that  they  lead  to  different  coarsening  kinetics. 

2.1.  Model  A 

In  this  model,  we  assume  that  the  YSZ  phase  has  a  very  small 
mobility,  and  consider  the  Ni-YSZ-pore  anode  as  a  system  of  three 
mobile  phases.  To  model  the  phase  evolution  in  a  three-phase  sys¬ 
tem  with  a  PFM,  three  OPs,  say  0i,  02,  and  03,  may  be  required. 
However,  by  choosing  the  OPs  to  be  the  volume  fraction  of  the 
corresponding  phases,  we  eliminate  one  of  the  OPs  through  the 
constraint  J^=10j  =  1.  The  governing  equations  then  reduce  to  two 
coupled  Cahn-Hilliard  equations  as 

fid).  RQf 

^  =V.M(01,02,03)V^,  (1) 

where  i  =  l,  2,  8F  is  the  free  energy  functional,  and  M  is  a  surface 
mobility  function,  described  below.  We  utilize  a  Ginzburg-Landau- 
type  free  energy  functional  [13]  of  the  form 

9  =  j  dV  f  Y  .  (2) 


excess  mobility  resulting  from  the  appearance  of  the  foreign  phase 
discussed  above.  Here,  we  set  Q  =  0.05  as  the  lower  cutoff  OP  value, 
and  Ch  =  1  -  Q  as  the  upper  cutoff  resulting  from  the  complemen¬ 
tary  value  of  the  lower  cutoff. 

To  quantitatively  link  the  simulation  predictions  to  exper¬ 
iments,  we  performed  an  asymptotic  analysis.  Following  Refs. 
[15,16],  we  obtain 


fcfiTN,L4 
8s  YsDs 


(4) 


where  kB  and  T  are  the  Boltzmann  constant  and  temperature.  Nv,  8S , 
ys,  and  Ds  are  the  atomic  number  density,  surface  diffusion  depth, 
surface  energy,  and  surface  diffusivity,  respectively.  L  and  r  are 
the  characteristic  length  and  characteristic  time  scale,  respectively. 
I<  is  approximately  unity  (equals  1.0147  for  the  simulations  pre¬ 
sented),  determined  by  analytical  integration,  and  it  compensates 
for  the  aforementioned  cutoff  in  the  mobility  function.  The  detailed 
derivation  will  be  published  in  a  forthcoming  paper. 


2.2.  Model  B 


In  this  alternative  model,  we  assume  YSZ  is  stationary  and  there¬ 
fore  can  be  treated  as  the  geometry  within  which  Ni  and  pore 
phases  evolve.  We  also  assume  that  the  Ni  phase  evolves  via  dif¬ 
fusion  along  Ni-pore  interfaces,  and  the  triple  junction  satisfies 
Young’s  equation  for  a  contact  angle  with  a  locally  flat  surface. 
The  dynamics  of  this  system  can  thus  be  described  by  a  single  OP 
Cahn-Hilliard  equation  with  two  complementary  boundary  condi¬ 
tions  (BCs):  the  contact-angle  BC  at  triple  junctions  and  no-flux  BC 
at  YSZ  interfaces,  where  the  OP  distinguishes  Ni  and  pore  phases. 
We  use  the  framework  of  the  Smoothed  Boundary  Method  (SBM) 
[19,20]  and  treat  the  contact-angle  BC  following  Warren  et  al. 
[17,18]. 

For  the  standard  Cahn-Hilliard  dynamics,  the  evolution  equa¬ 
tion  is 

^  =  V-MV^.  (5) 

The  chemical  potential  p,  can  be  calculated  as 
/x  =  <$3F/<50  =  3//3  0  -  s2  V20,  where  we  use  0  =  1  as  the  Ni  phase. 
The  free  energy  functional  in  Model  B  has  the  standard  form  that 
consists  of  a  double  well  potential  and  gradient  energy  penalty, 


/(0)  + 


(6) 


This  functional  allows  the  control  of  interfacial  energies  and 
thicknesses  of  each  type  of  interfaces  through  the  selection  of 
parameters  /q  and  ofy.  In  addition,  for  unequal  surface  tensions,  this 
functional  leads  to  more  computationally  efficient  evolution  equa¬ 
tions  as  compared  to  an  alternative  form  of  the  free  energy  in  Ref. 
[14].  The  variational  derivative  of  this  functional  can  be  calculated 
by  the  method  of  Lagrange  Multipliers  [14].  The  interfacial  energy 
Yij  =  ofy  y/(/q  +  kj)/3  and  interfacial  width  8y  =  /  y/(/q  +  /q)  for  a 

planar  interface  between  phases  i  and  j  can  be  found  from  the  equi¬ 
librium  solution.  The  model  permits  a  small  amount  of  the  third 
phase  to  appear  in  a  two-phase  interface  (as  discussed  in  Ref.  [14]), 
but  it  is  limited  to  less  than  3%  in  fraction. 

The  surface  mobility  function  is  taken  to  be 

3 

M(01,02,03)  =  Y  O) 

ij=i,i>j  qch  qch 


where /(0)  =  W02(1  -  0)2/4.  The  SBM  introduces  a  domain  param¬ 
eter  0,  which  in  this  case  distinguishes  regions  having  YSZ  (0  =  0) 
and  other  phases  (0  =  1)  as  well  as  the  YSZ  interfaces  (O<0<1). 
Like  an  ordinary  OP  in  a  phase-field  approach,  0  varies  continu¬ 
ously  across  the  interface;  thus,  the  unit  interface  normal  n  is  given 
by  n  =  V0/  |  V0|.  The  contact  angle  0  at  triple  junctions  dictated 
by  Young’s  equation  can  be  formulated  as  n  •  V0/  |  V0 1  =  -  cos  0. 

The  mechanical  equilibrium  at  the  triple  junction  corresponds 
to  an  extremum  of  the  free  energy,  or  8cP  =  0.  We  can  use  the 
planar  solution  of  the  volume  integration  part  of  this  extremum 
condition  to  find  a  useful  equality,  |  V0 1  =  y^2 f/e,  which  can  be 
substituted  into  the  aforementioned  Young’s  condition  to  derive 
the  SBM  contact-angle  BC  as  [17,21], 


V0  •  V0  =  -  |  V0|  cos  0 


(7) 


where  we  introduced  a  boxcar  function  YiqcJ^'  which  takes  a  This  contact-angle  BC  only  supplies  energy  on  interfaces  near  triple 
value  of  1  when  0/  is  between  Q  and  CB  and  0  otherwise,  to  avoid  junctions.  The  final  evolution  equation  derived  from  SBM  with  no- 
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Fig.  1.  A  comparison  of  equilibrium  contact  angles  at  triple  junctions  with  three  different  models;  phases  1,  2,  and  3  are  represented  in  green,  gray,  and  blue,  respectively: 
(A)  the  result  of  Model  A;  (B)  the  result  of  Model  B;  (C)  the  result  from  evolving  one  order  parameter  in  Model  A.  Note  that  the  triple  junction  no  longer  exists  at  equilibrium 
in  (C). 


flux  and  contact-angle  BCs  is 


W) 

dt 


=  v- 


•  (i/r V0)+ 


£ 


The  mobility  fuction  in  Model  B  is  given  by 

M(<p,  i /f)  =  MNi_pore  JJ(0X02(1  -  4>2))gW’ 

Q  0, 


(8) 


O) 


where  g(^)  =  ^6(10^2  - 15^  +  6)  is  introduced  to  control  the 
mobility  at  and  near  triple  junctions.  The  one-sided  interpolation 
function  g{\fr)  transits  smoothly  from  1  to  the  order  of  0.01  as  the 
domain  parameter  varies  from  1  to  0.5.  In  this  model,  the  mobility 
near  a  triple  junction  decreases  by  a  factor  of  10-6  as  ^  varies  from 
1  to  0.1  along  the  Ni-pore  interface. 


3.  Results  and  discussions 

3.1.  Contact  angles  at  triple  junctions 

The  mechanical  equilibrium  at  a  triple  junction  is  examined 
in  two  dimensions  with  713:723*712  =  1.9: 1 .4:2.2  to  validate  our 
models.  As  shown  in  Fig.  1,  the  contact  angles  in  both  Model  A 
and  Model  B  are  found  to  obey  the  Young’s  equation  at  equi¬ 
librium  with  different  assumptions:  in  Model  A,  YSZ  surface 


deforms  near  the  TPB  while  in  Model  B,  the  YSZ  is  not  allowed  to 
deform. 

A  model  one  may  consider  is  Model  A  (with  two  OPs)  where 
one  of  the  OPs  is  held  constant  after  proper  initialization.  Such  a 
model  has  been  proposed  and  applied  by  others  [8,1 3  ].  We  find  that 
this  model  would  only  work  when  the  mobile  phase  completely 
wets  the  stationary  phase  as  in  Ref.  [13],  but  not  when  three-phase 
boundaries  are  involved.  In  a  PFM,  the  interfacial  energy  is  con¬ 
trolled  by  the  bulk  and  gradient  energy  coefficients,  and  an  interface 
possesses  the  prescribed  interfacial  energy  when  the  OP  interfacial 
profile  is  in  equilibrium  (or  near  equilibrium  when  the  interface 
is  evolving).  In  Eq.  (1),  if  one  OP  is  held  constant  as  proposed  in 
[8],  other  OPs  cannot  develop  their  equilibrium  profiles.  This  leads 
to  contact  angles  at  triple  junctions  that  are  inconsistent  with  the 
solutions  of  Young’s  equation.  The  full  analysis  of  this  problem  is 
beyond  the  scope  of  this  paper.  Here,  we  simply  demonstrate  that 
the  Young’s  condition  is  not  met  in  Model  A  with  one  OP  held 
constant,  as  shown  in  Fig.  1(C). 

3.2.  Effects  of  coarsening  on  micro  structure  and  TPB  density 

Microstructures  of  SOFC  anodes  can  be  experimentally  recon¬ 
structed  [5].  We  use  an  experimentally  obtained  Ni-YSZ  anode 
microstructure  of  dimension  3.48  by  3.48  by  3.48  p,m3  (Fig.  2, 
left)  as  the  initial  condition  for  our  coarsening  simulations.  For 
convenience,  hereafter  subscripts  1,  2,  and  3  are  used  to  repre¬ 
sent  the  Ni,  YSZ  and  pore  phases,  respectively.  Parameters  My 
and  yij  are  determined  based  on  the  physical  properties  of  the 
Ni-YSZ  system.  Surface  tensions  of  Ni  and  YSZ  are  determined  to 


Fig.  2.  A  comparison  of  the  three  phase  anode  before  (left)  and  after  (right)  coarsening  for  100  h  for  the  reference  case  using  Model  B.  The  Ni,  pore,  and  YSZ  phases  are 
represented  in  green,  blue,  and  semi-transparent,  respectively  (with  volume  fractions  24.2%,  20%,  and  55.8%,  respectively).  The  TPB  density  reduced  by  23%  after  coarsening. 
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Fig.  3.  A  comparison  of  the  TPB  reduction  between  Model  A  and  Model  B  for  the 
reference  case,  as  well  as  for  three  different  Ni-YSZ  interfacial  energies  using  Model 
B. 


be  y i3  =  1.9  J  m-2  and  y2 3  =  1  .4J  m- 2  [22],  and  surface  diffusivity  of 
them  are  D13  =3  x  10-10m2s-1  [10,11]  and  D23  =3  x  10_16m2s_1 
[12].  We  estimate  the  interfacial  tension  of  Ni-YSZ  interface  to 
be  yi2  =  1.5Jm-2  [23,24]  and  use  Du  =  3  x  10-20  m2  s_1,  respec¬ 
tively.  Therefore,  we  choose  Mi3:M23:Mi2  =  1 : 1 0— 6 : 1 0— 10  and 
Yi3:Y23-Yi2  =  1.9:1.4:1.5  as  the  reference  case,  using  a  mobility 
scale  of  3  x  1 O-10  m2  s_1  and  an  energy  scale  of  1 J  m-2. 

A  comparison  of  the  anode  microstructure  before  and  after 
coarsening  for  1 00  h  using  Model  B  is  shown  in  Fig.  2.  The  evolution 
of  the  microstructure  appears  to  be  confined  within  a  short  range. 
Note  there  are  a  few  regions  where  YSZ  structure  appears  to  change, 
but  they  are  visualization  artifacts.  As  shown  in  Fig.  3,  Model  B  pre¬ 
dicts  that  the  three-phase  boundary  (TPB)  length  becomes  nearly 
constant  after  coarsening  for  tens  of  hours,  while  Model  A  indi¬ 
cates  that  the  TPB  continues  to  decrease  slowly  after  a  rapid  early 
evolution.  These  results  suggest  that  Ni  phase  can  be  trapped  in 
metastable  states  if  YSZ  is  immobile.  Interestingly,  the  TPB  reduc¬ 
tion  during  the  early  stage  of  coarsening  are  very  similar  for  both 
models. 

For  the  reference  case,  the  overall  TPB  reduction  after  100- 
h  coarsening  was  found  to  be  23%  for  Model  B.  A  preliminary 
experimental  result  for  the  corresponding  case  was  found  to  be 
approximately  26%  [25],  which  is  within  the  uncertainties  asso¬ 
ciated  with  the  simulation.  Sources  of  uncertainties  include  the 
values  of  diffusivities  (which  change  the  time  scale),  the  interfa¬ 
cial/surface  energies  (which  change  the  morphological  evolution 
as  well  as  the  time  scale),  statistical  variation  of  the  microstructure 
in  the  selected  volume,  and  possible  initial  smoothing  of  very  small 
unresolved  structures.  TPBs  in  Model  A  reduced  by  44%  after  33-h 
coarsening,  and  thus  it  overestimates  the  reduction  significantly. 
This  overestimation  is  mainly  due  to  the  fact  that  the  mobility  func¬ 
tion  we  used  can  introduce  excess  mobility  at  triple  junctions  even 
with  the  cutoff,  especially  when  the  mobilities  differ  by  orders  of 
magnitude  as  in  the  present  case. 

3.3.  Effects  of  contact  angles  at  triple  junctions  on  anode  stability 

A  parametric  study  of  three  different  Ni-YSZ  interfacial  ener¬ 
gies  (leading  to  three  different  sets  of  contact  angles)  is  performed 
using  Model  B.  As  shown  in  Fig.  3,  the  general  trend  of  the 
coarsening  behavior  remains  the  same  for  all  parameter  sets. 
However,  TPBs  stabilize  at  an  earlier  time  with  less  reduction 
when  the  Ni-YSZ  interfacial  energy  is  smaller,  which  corre¬ 
sponds  to  a  smaller  contact  angle  of  Ni  on  YSZ.  Therefore,  such 
a  case  is  found  to  be  advantageous  in  preserving  TPBs  in  the 
Ni-YSZ  anode.  This  result  indicates  that  there  is  an  opportu¬ 


Fig.  4.  The  evolution  of  the  tortuosity  factor  and  the  interfacial  area  reduction  of 
the  pore  phase  predicted  by  Model  B. 


nity  for  improving  the  long-time  performance  of  SOFC  anodes 
by  controlling  the  surface/interfacial  energies  of  materials,  which 
can  be  accomplished,  for  example,  by  doping  Ni  with  additives 
[22]. 

3.4.  Effects  of  coarsening  on  pore  phase  tortuosity 

The  pore  phase  in  SOFC  anodes  facilitates  the  transport  of 
hydrogen.  Tortuosity  through  the  pore  microstructure  allows  us 
to  characterize  the  overall  transport  property,  and  it  can  be  used  in 
modeling  the  electrochemical  performance  of  anodes  via  macro- 
homogeneous  models  [26].  However,  tortuosity  is  difficult  to 
determine  experimentally.  Therefore,  using  the  microstructures 
evolved  by  Model  B,  we  examine  the  evolution  of  the  tortuos¬ 
ity  factor  determined  as  described  in  [5,27],  averaged  over  three 
orthogonal  directions,  and  the  interfacial  area  of  the  pore  phase 
over  lOOh  of  coarsening.  As  can  be  seen  in  Fig.  4,  the  tortuosity 
variation  is  correlated  to  the  interfacial  area  changes  during  the 
late  evolution;  however,  this  is  not  necessarily  true  in  the  early 
stage.  The  transport  in  the  pore  phase  is  thus  not  solely  affected  by 
the  average  length  scale,  which  is  defined  here  as  the  inverse  of  the 
surface  area  per  unit  volume.  Our  results  suggest  that  coarsening 
for  a  significant  amount  of  time  may  have  positive  effects  on  the 
gas  transport  in  the  pore  phase. 

4.  Conclusions 

In  conclusion,  we  present  the  framework  of  a  phase-field 
approach  to  quantitatively  model  Ni  coarsening  in  SOFC  Ni-YSZ 
anodes.  Two  models  are  presented— Model  A  permits  limited  YSZ 
mobility  while  Model  B  assumes  YSZ  to  be  immobile.  Although  our 
models  contain  approximations  and  the  simulation  results  may  be 
affected  by  uncertainty  in  material  properties,  we  find  a  reason¬ 
able  agreement  of  the  TPB  reduction  trend  between  coarsening 
experiments  and  simulations  using  Model  B.  Model  A,  however, 
overestimates  the  TPB  reduction  due  to  the  evolution  of  the  YSZ 
structure  induced  by  an  excess  mobility  at  triple  junctions  con¬ 
tained  in  the  model.  Our  simulations  show  that  the  major  portion 
of  TPB  reduction  occurs  in  the  early  stage  of  coarsening,  and 
stability  of  TPBs  is  observed  provided  that  YSZ  is  stationary.  Fur¬ 
thermore,  we  demonstrate  that  smaller  Ni-YSZ  interfacial  energy 
can  enhance  the  stability  of  the  TPB  density.  These  models  can  be 
extended  to  simulate  other  operating  conditions  of  SOFCs.  More¬ 
over,  the  framework  of  Model  B  can  be  applied  to  solve  other 
physical  problems  when  the  wetting  condition  plays  an  important 
role. 
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